import numpy as np
import random
import time


def hecheng2(time, w, fai, Aw, bl):  # 由输入参数合成地震波  旧方法
    nt = len(time)
    npl = len(w)
    atj = np.ones((len(w), len(time)))
    for i in range(npl):  # 频率点
        inner_fai = fai[i]
        inner_w = w[i]
        inner_Aw = Aw[i]
        for j in range(nt):  # 时间点  某频率的相位固定
            atj[i, j] = np.sin((inner_w * (time[j]) + inner_fai)) * inner_Aw * bl[j]
    at = np.sum(atj, axis=0)  # 各个频率在各个时间点的分量叠加

    return at


def hecheng(time, w, fai, Aw, bl):  # 由输入参数合成地震波   想要快速 就要能够对w Aw 数组形式进行计算
    """
    合成地震波
    :param time: 时间序列
    :param w:    频率
    :param fai:  相位
    :param Aw:   幅值
    :param bl:   外包络
    :return:
    """
    # 纯数组操作
    (w_xlen, w_ylen) = w.shape if w.ndim == 2 else (len(w), 0)
    (Aw_xlen, Aw_ylen) = Aw.shape if Aw.ndim == 2 else (len(Aw), 0)
    nw = len(w)
    nt = len(time)
    if w_ylen:  # 如果w是多维数组：
        # print('w是多维数组')
        "=======================将原始维度拓展到三维======================="
        ww = w.reshape(w_xlen, 1, w_ylen)
        tt = np.dot(np.ones((w_xlen, 1)), time.reshape(1, nt)).reshape(w_xlen, nt, 1)  # 拓展到二维
        faiffai = np.dot(fai.reshape(len(fai), 1), np.ones((1, nt))).reshape(w_xlen, nt, 1)  # 拓展到二维
        blbl = np.dot(np.ones((w_xlen, 1)), bl.reshape(1, nt)).reshape(w_xlen, nt, 1)  # 拓展到二维
        AwAw = np.dot(Aw.reshape(len(Aw), 1), np.ones((1, nt))).reshape(w_xlen, nt, 1)  # 拓展到二维
        "====================拓展维度后，进行计算 ===================="
        atj = ww * tt + faiffai
        atj = np.sin(atj)  # 求sin
        atj = atj * blbl * AwAw
        at = np.sum(atj, axis=0)  # 各个频率在各个时间点的分量叠加 sum到二维
    elif Aw_ylen:  # 如果Aw是多维数组：
        # print('Aw是多维数组')
        "=======================将原始维度拓展到三维======================="
        ww = np.dot(w.reshape(nw, 1), np.ones((1, nt))).reshape(nw, nt, 1)  # 拓展到二维
        tt = np.dot(np.ones((nw, 1)), time.reshape(1, nt)).reshape(nw, nt, 1)  # 拓展到二维
        faiffai = np.dot(fai.reshape(nw, 1), np.ones((1, nt))).reshape(nw, nt, 1)  # 拓展到二维
        blbl = np.dot(np.ones((nw, 1)), bl.reshape(1, nt)).reshape(nw, nt, 1)  # 拓展到二维
        AwAw = Aw.reshape(nw, 1, Aw_ylen)
        "====================拓展维度后，进行计算 ===================="
        atj = ww * tt + faiffai
        atj = np.sin(atj)  # 求sin

        atj = atj * blbl * AwAw   # 三维乘法
        # 对于维数一样的张量相乘。A.shape =（b, m, n), B.shape = (b, n, k).
        # numpy.matmul(A, B)
        # 结果shape为(b, m, k)
        # 这里要求第一维度相同。后两个维度能满足矩阵相乘条件。

        at = np.sum(atj, axis=0)  # 各个频率在各个时间点的分量叠加 sum到二维
    else:  # 正常数组
        # print('正常数组')
        atj = w.reshape(nw, 1) * time + fai.reshape(len(fai), 1)  # 乘上w和t
        atj = np.sin(atj)  # 求sin
        atj = atj * Aw.reshape(len(Aw), 1) * bl  # 乘上包络和幅值系数
        at = np.sum(atj, axis=0)  # 各个频率在各个时间点的分量叠加
    return at


if __name__ == '__main__':
    t = np.arange(0, 20, 0.02)  # 人造地震动时间序列
    dltf = 0.05
    w = np.arange(0.15 * 2 * np.pi, (25 + dltf) * 2 * np.pi, dltf * 2 * np.pi)  # 角频率  %频率变化范围0.15Hz~25Hz
    fai = np.zeros(len(w))  # 随机相位序列
    n = 5
    Aw = np.ones(shape=(len(w), n))
    bl = np.random.randint(1, 9, len(t))
    s = time.time()
    for i in range(n):
        hecheng2(t, w, fai, Aw[:,i], bl)
    e = time.time()
    t1 = e-s
    print(t1)
    s2 = time.time()
    b = hecheng(t, w, fai, Aw, bl)
    e2 = time.time()
    t2 = e2-s2
    print(t2)
    # print(a == b[:,1])
    print(t1/t2)
    print(t.shape)
    print(w.shape)
    print(Aw.shape)